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Continuously sell-similar (CSS) solutions for the gravitational collapse of a spherically symmetric 
perfect fluid, with the equation of state p — Kp, with < k < 1 a constant, are constructed 
numerically and their linear perturbations, both spherical and nonspherical, are investigated. The 
/ = 1 axial perturbations admit an analytical treatment. All others are studied numerically. For 
intermediate equations of state, with 1/9 < k < 0.49, the CSS solution has one spherical growing 
mode, but no nonspherical growing modes. That suggests that it is a critical solution even in 
(slightly) nonspherical collapse. For this range of k we predict the critical exponent for the black 
hole angular momentum to be 5(1 + 3k)/3(1 + k) times the critical exponent for the black hole mass. 
For « = 1/3 this gives an angular momentum critical exponent of /j, ~ 0.898, correcting a previous 
result. For stiff equations of state, 0.49 J$ k < 1, the CSS solution has one spherical and several 
nonspherical growing modes. For soft equations of state, < k < 1/9, the CSS solution has 1+3 
growing modes: a spherical one, and an I — 1 axial mode (with m = —1, 0, 1). 
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I. INTRODUCTION 

An isolated system in general relativity ends up in one of three stable final states: a black hole, a star, or complete 
dispersion. The phase space of isolated systems in general relativity is therefore divided into basins of attraction: 
each initial data set must end up in one of the stable end states. The study of the boundaries between the basins of 
attraction, in particular of the boundary between the black hole and dispersion end states, began with the pioneering 
work of Choptuik Jlj] , and is now an active field in classical general relativity. 

Initial data near the black hole/dispersion threshold evolve through a universal intermediate state before dispersing 
or forming a black hole. This intermediate attractor has higher symmetry, as a spacetime, than the generic solution. 
Often it is self-similar. Close to the threshold, but on the collapse side, the mass of the final black hole then scales as 
a universal power of the distance of the initial data to the black hole threshold. Universality, self-similarity and the 
critical exponent for the black hole mass have given rise to the name (type II) "critical phenomena in gravitational 
collapse" . For a review see 0] . 

Critical phenomena of this type are explained by the existence of a solution that is self-similar, regular, and has 
exactly one growing perturbation mode, such that for one sign of the growing mode the solution veers towards black 
hole formation, and for the other towards collapse. Such a solution is called a (type II) critical solution. From 
a dynamical systems point of view, a critical solution is an attractor within the black hole threshold, which is a 
hypersurface of codimension one. Within the complete phase space, it is therefore an attractor of codimension one. 
All solutions that start near the black hole threshold, but not necessarily near the critical solution itself, are funneled 
through this intermediate attractor. This funneling process explains both universality and the self-similar nature of 
the intermediate attractor explains scaling. The critical exponent in the power-law scaling of the black hole mass can 
be shown to be the inverse of the Lyapunov exponent of the critical solution's one growing perturbation mode |3j. 

Here we concentrate on one class of matter models coupled to general relativity, perfect fluids with the linear 
barotropic equation of state p — up, where p is the pressure, p is the total energy density measured in the rest frame, 
and k is a constant in the range < k < 1. The spherically symmetric fluid with k = 1/3, corresponding to an 
ultrarelativistic gas, was one of the first systems in which critical phenomena were found 0] . These results were later 
extended to the range < k < 1 p-R]. We now ask what happens when we allow small deviations from spherical 
symmetry. 

For a sample of values of k in the range < k < 1, we construct a regular, continuously self-similar (CSS), and 
spherically symmetric solution, and then investigate its linear perturbations to see how many growing perturbation 
modes it has. It is already known, from both perturbative [g,|6| and non-perturbative J7|] calculations that these 
solutions have exactly one growing mode among their spherically symmetric perturbations, which makes them critical 
solutions in spherical symmetry. Here we examine their nonsphcrical perturbations. In a previous Rapid Communi- 
cation @] we examined the particular case k — 1/3. Here we generalize this investigation to all values of k in the 
range < k < 1. We also describe our numerical methods and results in much more detail. Finally, we correct two 
incorrect assumptions in H , namely that the I = 1 axial perturbations obey the same type of equation as the / > 2 
ones, and that the I = 1 polar perturbations are pure gauge. As it happens, correcting these errors does not affect the 
overall conclusion of Ref. 0, namely that the critical solution for k = 1/3 has no growing nonspherical perturbation 
modes. 

As the background is spherically symmetric, the perturbations can be separated into spherical harmonics labelled 
by / and m, and can be split further into axial and polar parts. The perturbation equations we use are those derived 
in U , restricted here to a CSS background and the particular equation of state p = up. Surprisingly, the nonspherical 
perturbation equations are much harder to solve numerically than in the equivalent problem for the massless scalar 
field fllCf . One difficulty is that we are dealing with a one-parameter family of background solutions, whose limits 
k = and k = 1 are not regular members of the family. The other difficulty is that in the I > 2 perturbations both 
light cones and sound cones play a dynamical role, while they coincide for the scalar field. This gives rise to weak 
solutions (in the sense of hyperbolic equations) that we need to discard. Because of these problems, our final choice 
of numerical approach is to discretize the perturbation evolution equations in space but not in time. We then look 
directly for eigenvectors and eigenvalues (modes) of the finite difference equations. 



The plan of the paper is this: in section y we discuss the general perturbation framework of Ref. ||, and the CSS 
background solutions. Section [II discusses the I = 1 axial perturbations. Their spectrum can be calculated in closed 



form. Based on this, we correct the value of the angular momentum critical exponent stated in Ref. |11|. All other 



perturbations require a numerical treatment and are discussed in Section [V . Details of the numerical difficulties and 
numerical methods, however, are given in the appendixes. Section [V] summarizes our results. 

II. BACKGROUND SOLUTION 

A. Perturbations of a spherically symmetric perfect fluid 

We shall examine the linear perturbations of a spherically symmetric and continuously self-similar (CSS) perfect fluid 
spacetime. For this purpose we use the restriction to a CSS background of a general framework for the perturbations 
of time-dependent spherically symmetric perfect fluid solutions that was presented in [jpj. In this formalism, the 
spacetime manifold is written as the product M = M 2 x S 2 , where S 2 is the 2-sphere and M 2 is a 2-dimensional 
manifold, the "rt-plane", with a boundary r — 0. The coordinates in in M 2 are denoted by x A , and the coordinates in 
S 2 by x a . Coordinates in all of M are collectively denoted by x M . The general spherically symmetric metric becomes, 
in this notation, 

g M „ = diag (gAB,r 2 j ab ) , (1) 

where r 2 is a scalar function on M 2 , and 7 a & is the unit metric on S 2 . The spherically symmetric perfect fluid 
stress-energy tensor can be written in the same notation as 

T^ = diag (p u A u B +pn A n B ,p r 2 j ab ) , (2) 

where p and p are the (total) energy density and pressure in the fluid rest frame, u M = (u , 0) is the fluid 4- velocity, 
and n A is the outward pointing unit vector normal to u A in M 2 . The held equations in this framework are covariant 
in M 2 . The two extra dimensions in S 2 appear through the scalar r. 

The perturbations of a spherically symmetric background decompose naturally into polar and axial parity, and into 
spherical harmonic angular dependencies, for different I and m. The equations of motions are the same for all values 
of — I < m < I, for given I > 0. The cases I = 0, 1 = 1 (polar and axial), and I > 2 (polar and axial) are all qualitatively 
different and need to be treated separately. Tensors in M, including the perturbations, are written as products of 
tensors in M 2 with tensors in S 2 . All necessary tensors in S 2 are built from the scalar spherical harmonics Y/ m on 
S 2 , their covariant (with respect to jab) derivatives, and the covariantly constant antisymmetric tensor e ao . The final 
equations for the perturbations are again covariant equations for tensors on M 2 . Their angular dependence comes in 
through terms such asl/r and 1(1 + l)/r 2 . 

In the next step, linear combinations of the perturbations are found that do not change to linear order under 
inhnitcsimal coordinate transformations, in cither M 2 or S 2 . The perturbation equations can be rewritten in terms 
of these gauge-invariant perturbations alone. 

In a further step, all perturbation tensors in M 2 are split into frame components with respect to the orthonormal 
frame (u A , n A ). This "scalarization" replaces covariant derivatives of tensors with partial derivatives of scalars. These 
derivatives are also decomposed into their frame / = u f,A and /' = n f,A- The perturbation equations are now 
scalar equations written without reference to a particular coordinate system. In this sense they are covariant, as well 
as linearly gauge- invariant. Note that the frame derivatives are not partial derivatives, and do not commute. The 
advantage of using these derivatives is that / ± c s f are derivatives along radial matter characteristics, and / ± /' 
derivatives along the radial light rays. On the other hand, some constraint-type perturbation equations for I = and 
I = 1 are most naturally written using the derivative D = rdf/dr along polar slices. 

In a final step, we rescale the perturbation variables by ^-dependent powers of r so that they are either O(l) at the 
origin and even functions of r, or else 0(r) and odd. The equations are brought into first-order form by treating first 
derivatives such as / and /' as independent variables where necessary. 

In the remainder of this section, we introduce coordinates that are adapted to a CSS background, and review how 
the background solution is defined, and constructed numerically. 



B. Continuously self-similar background solution 

Although our perturbation variables are linearly gauge-invariant, we have to coordinate system on the background. 
A standard choice of coordinates x A is to use r as a coordinate (radial gauge), and make the second coordinate t 
orthogonal to it (polar slicing) . Then qab takes the form 

g tt = -a 2 , g rr = a 2 , g rt = 0, (3) 

with a and a functions of r and t. Radial light rays are governed by the combination g = a/a. There is a remaining 
gauge freedom t — > t'(t), and we fix it by setting a — 1 at r = for all t. 

Based on radial-polar coordinates, we now introduce coordinates that are adapted to self-similarity, while retaining 
polar slicing. We define new coordinates x and r by 

r = sxe~ T , t = -e~ T , (4) 

with s > a constant. We have assumed that t < 0, and have chosen signs so that r increases as t increases. Note 
that r — > oo as t — > 0_. Partial derivatives transform as 

f,t = e T (/,r + xf lX ), f, r = s- x e T f, x . (5) 

The metric g&B in these coordinates becomes 

— 2t/ 2 i 2 2 2\ — 2t 2 2 — 2t 2 2 / '^\ 

g rr = e (—a + s x a ), (j^ = e s a , <? r:E = — e s xa . (6) 

The spacetime is continuously self-similar (with homothctic vector —d/dr) if and only if a and a depend only on 
x but not on r. r has two different interpretations. On the one hand, it is a time coordinate in the sense that its 
level surfaces are spacelike. But — r is also the logarithm of spacetime scale, in the sense that proper distances are 
proportional to intervals Ax and At times a factor of e~ r . In a self-similar spacetime, larger r therefore means 
structure on a smaller scale. The point r = 0, t = 0, or r = oo, is by construction a curvature singularity, unless the 
spacetime is flat. 

In the solutions we consider here the matter is a perfect fluid with density p and pressure p = up, with k a constant. 
This equation of state is the only one compatible with exact self-similarity. We impose CSS in the metric by making 
the ansatz a = a(x) and a = a(x). Imposing self-similarity on the spacetime, we find from the Einstein equations 
that Airp = e 2l ~p{x) and V = u r^/n B r.B = ^(a;). Here n = — e bu b is the outward-pointing unit spacelike 
vector normal to u . The constant s that was introduced above is now chosen so that the surface x — 1 is a matter 
characteristic. It is then the past sound cone of the singularity. 

The background solution that we want to use is completely defined by the assumptions of i) continuous self-similarity, 
ii) spherical symmetry, iii) analyticity at the center x = 0, and iv) analyticity at the past sound cone x = 1. The 
background equations resulting from the CSS ansatz are given in Appendix |Aj. The CSS ansatz reduces the two 
Einstein equations and two matter equations that are needed in spherical symmetry to one algebraic equation for a 
three ordinary differential equations in x for p, V and g = a/a. 

There are two boundary conditions for a CSS solution at x = 0. The gauge condition a(0) — 1 becomes g(0) = 1. 
From regularity of the matter velocity, CSS, and matter conservation one can derive that 

lim — = - —. 2 . . (7) 

x^o sgx 3(1 + k) 

Imposing this is the regularity condition iii) at the center. Note that this limit would also hold for a CSS fluid in a 
flat spacetime. 

One boundary condition at X = 1 is the gauge condition D(l) = 0, which makes x = 1 is the sound cone. This 
condition determines the value of the constant s. The regularity condition iv) at the sound cone is 

y/HSi(l) = S 2 {1), (8) 

which is the vanishing of the term in the equations that is divided by D (see Appendix |A|) . 

Once we have solved the boundary value problem in < x < 1, we can analytically continue the solution through 
x = 1 (which is a regular singular point of the equations) and then continue the solution by evolving the ODEs to 
larger x. (In numerical terms, analytic continuation is implemented by polynomial extrapolation.) We go to the light 
cone and a little beyond. The light cone is at a value of x that depends on k. The ODEs are regular there. 



Our numerical method for imposing analyticity at X = 1 and x — is to just impose the algebraic boundary 
conditions there, and to use centered differences everywhere else, without using an explicit power-law expansion 
about the singular points. The first (numerical) derivative of the fluid density and velocity profiles obtained with 
this method has a small discontinuity at x = 1 that first appears at k ~ 0.7 and increases with n. Results for 
k > 0.7 therefore have a source of numerical error over and above the one arising in the numerical evolution of the 
perturbations. 



III. AXIAL L = 1 PERTURBATIONS - ANALYTICAL TREATMENT 



A. Equation of motion 

In this section we discuss the axial I = 1 perturbations of the continuously self-similar perfect fluid critical solution. 
This leads us, from general arguments given in |11|, to a prediction for the scaling of black hole angular momentum 
in critical collapse. Note that this section differs from the rest of the paper in presenting purely analytical results. 

The axial I — 1 perturbations contain a single matter degree of freedom, and no gravitational waves. The gauge- 
invariant fluid velocity perturbation, 0, obeys the autonomous equation of motion 



(Pr 2 pu A ) lA = 0. 



(9) 



This is just a transport equation along the background fluid flow. All axial metric perturbations are encoded in 
a gauge- invariant scalar II. For I = 1, II is obtained from (3 by quadrature, and Eq. (ft) describes the dynamics 
completely. For I > 2, II obeys a wave equation with a source term proportional to j3. (In |11| , it was incorrectly 
assumed that this is true also for I = 1.) 

As we shall see now, the complete mode spectrum of (3 can be obtained analytically for all I. For I = 1, we then 
have the complete dynamics. For I > 2, the spectrum of is also known analytically, but that of the homogeneous II 
modes must be calculated numerically. Here we obtain the modes of for general I > 1, and restrict to I = 1 at the 
end. 

The perturbed fluid velocity is regular at r = if is 0(r l+1 ) there. We therefore define a rescaled variable 
(3 = r~'' +1 '/3 that is even in x and generically nonzero at x = 0. The perturbed stress-energy tensor remains self-similar 

o _ 

if depends on r as e~ T at constant x (for any I). We therefore define a rescaled variable (3 = e~ lT = e T (sx)~^ +1 " > 0. 
This final variable obeys the equation 



■A, - \ X +J- ) '• 



1 + K 



2 -f | X+— ) (lli/M., 



+ l / + r,(i + -£-)>, o 



(10) 



This is of the form 



T + xA(x){3 x + B{x)(3 = 0, 



(11) 



where A and B are regular, even, strictly positive functions of x, so that A(x) = Aq + A2X 2 + 0(x 4 ), and similarly 

o 

B(x) = Bq + B 2 x 2 + 0(x ). We look for solutions (3 that are regular even functions of x. 



B. Analytic calculation of the mode spectrum 



Using the method of characteristics, the general solution of (|11|) can be written as 



I3{x, t) = e- BoT exp 



x B(x) - B 
xA(x) 



dx F 



xexp 



A{x) - Aq 
xA(x) 



dx — Aqt 



(12) 



where F{z) is a free function that is determined by the initial data. Note that the two definite integrals exist and are 

o 

0(x 2 ). For regular, even initial data, with f3(x,0) = Fq + 0(x 2 ), we have F(z) = Fq + 0(z 2 ). 

o 

We can now read off that the late-time behavior of the solution f3(x, r) as r — * 00 is 



(3(x, t) = e- BoT [F Q + 0(x 2 ) + O(e~ 2A0T )] 



(13) 



as r — > oo. For generic regular initial data, Fq does not vanish, and the solution decays as e~ B ° T at late times. One 
might have expected that the growth exponent depends on details of the background, but in fact it depends only on 
the background at the center. The physical reason for this is that the transport equation transmits information only 

from the center outwards. For example, we can see from (|12|) that the solution J3 for initial data that vanish in a 
neighborhood of x — is strictly zero at any fixed x at sufficiently large r. 

Surprisingly again, we can evaluate Bq in closed form, even though the background solution away from the center 
is known only numerically. From matter conservation and the assumption of continuous self-similarity we have the 
regularity condition (J7|). When we also take into account that p is an even function of x and that v is an odd function, 
we find from ( JIOJ ) that 

3(1 + k) V ; 

This analytic result is in perfect agreement with the late time behavior of numerical evolutions of generic initial data 

for/3. 

We now show that the mode spectrum is discrete. To look for modes, we make the ansatz 

kx,r) = e Xr f(x) (15) 

and obtain 

x Tx = -^r ! (16) 

By expanding this equation in powers of x around x = and comparing coefficients, we see that if /(0) ^ 0, and / is 
to be a regular even function of x, we must have A = —Bq. This is precisely the mode that dominates the late-time 
behavior (|13|), and we may call it fo(x). One obtains it as a power series in a neighborhood of x — 0, and then by 



solution of the ODE (16). We normalize Jq by setting /o(0) = 1. We now subtract a suitable multiple of this mode 

o 

from the solution (3 to obtain something that is 0(x 2 ) at x = 0: 

0™ ( x> T ) = fifa T ) _ Foe -Bo r fo{x) (17) 

o o 

where j3(x,0) = F -\-O(x 2 ) as before. By construction /?( 2 )(a?,0) = 0(x 2 ). Expanding ( |l2| ) to the next order, we have 

^ 2 \x,T) = e-^' +2A ^ [F 2 +0(x 2 ) + 0(e~ 2A ° T )]x 2 (18) 



We then obtain a mode f%(x) by solving (16) with A = —(-Bo + 2A ). We normalize it as fz(x) — x + 0(x ). 
Continuing in this way, we can strip off a sequence of modes decaying with A = —(Bo + 2uAq) for n = 0, 1, 2, . . .. 

o 

These modes are 0(x 2n ) at the origin, and so form a complete basis for smooth functions j3. Therefore the entire 
spectrum is discrete. 

From (|l0|), using again (R]), we find 

1 + 3k ,_ . 

3(1 + K) 

We can now write down the complete spectrum for all values of k and I, labelled by the index n — 0, 1, 2, It is 

w , n Df a o,^ 2(l-3«)-(l + 3/s)(i + 2n) ,„. 

X(k, I, n) = -B (k, I) - 2uA (k) = -i ) '-. 20 

3(1 + K) 

From this formula we can read off that alH > 2 modes decay for all k in the range < n < 1 . All I = 1 modes also 
decay for k > 1/9, but for k < 1/9 there is exactly one growing I = 1 mode (the n — mode). A for the dominant 
(n = 0) I = 1 mode is relevant for angular momentum scaling. It is 

1 - 9k 

Al = SM.L V 21 

3(1 + K) 

The dominant (n — 0) mode is plotted for I = 1 ... 5 in Fig. nL 




0.2 0.4 0.6 0.8 1 

o 

FIG. 1. The J3 equation describes axial fluid velocity perturbations. The plot shows the growth exponent A, which is real, 
against k. We plot X(k, l,n) given in Eq. fed). The thick line is for I — 1 and n = 0, where (3 is the only perturbation. (It is 
positive for k < 1/9.) Below, from top to bottom I — 2 . . . 5, for the leading mode, with n = 0. 



C. Angular momentum scaling 



In pT| we derived a general formula for the critical exponent /i governing black hole angular momentum in critical 
collapse, namely 



2- Ax 



(2-Ai) 7 . 



(22) 



(This formula corrects a misprint in Equation (11) of Reference [111.) Here Ao is the Lyapunov exponent for the 
spherical mode, which is real and positive, and Ai = A(k, 1,0) is the Lyapunov exponent for the I = 1 axial per- 
turbations, which is real. 7 = I/Aq is the critical exponent for the black hole mass. The derivation of this formula 
assumes that the critical solution has precisely one growing perturbation mode, which is spherically symmetric, while 
all nonspherical perturbation modes decay. The analytical result in this section and the numerical results in the 
following sections show that this assumption holds in the r ang e 1/9 < k < 0.49. 

For this range of equations of state, putting the results (|lj) and ( f22| ) together, we obtain the analytical prediction 



h(k) = 



5(1 + 3k) 
3(1 + k) 



7(«), 



< k < 0.49, 



(23) 



for the angular momentum exponent [A, given the mass exponent 7. 

This result corrects the value of \i given in |11| for the case k = 1/3, based on an incorrect value of Ai. The correct 
value Ai = —1/2 for k — 1/3 gives a real critical exponent /x = 57/2 ~ 0.898 for the black hole angular momentum. 
Note that the incorrect value given in pTj was complex, which was expected to give rise to oscillations in the direction 
of the black hole angular momentum as the black hole threshold is approached. However, the correct value of Ai is 
real, and so the angular momentum scaling is a power law, like the mass scaling. 



IV. ALL OTHER PERTURBATIONS - NUMERICAL TREATMENT 



A. General aspects of the perturbation equations 



Before we discuss the perturbation equations of motion in detail, it is useful to discuss the general form of the 
perturbations of a spherically symmetric and CSS solution. As we have seen, we can choose dependent and independent 
variables for the background so that the background is given by a number of functions Z(x,t) — Z* (x) of a single 
self-similarity variable x. It follows that the equations of motion of the linear perturbations of this background, when 
written in first-order form in suitable variables, are of the general form 



u, T = A(x)u, x + B(x)u, (24) 

u, x = C{x)u. (25) 

Here u(x,t) is a vector of perturbations, and A, B, C are matrices that depend on the background. We shall call 
the coupled partial differential equations that contain r-derivatives evolution equations, and the coupled ordinary 
differential equations in x constraints. Not all variables u need obey both types of equation - some variables u obey 
only an evolution equation, others only a constraint equation, and yet others both. For most of the perturbations, 
we shall be able to use a free evolution scheme in which the only constraints are trivial ones of the form u\ jX = 112 
introduced by writing a wave equation in first-order form. For the spherical perturbations and the polar I = 1 
perturbations, however, it is unavoidable to solve nontrivial constraint equations. 
The perturbation equations ( |24| , |25| ) admit solutions of the form 

u(x,t) = e XT u x {x). (26) 

Both A and u\{x) are in general complex, but because the coefficients A, B and C are real, the solutions form complex 
conjugate pairs. With A = A ± iui, the general real solution is 

Re[Ce iS u(x, r)] = Ce Ar [cos{lut + S) Reu x {x) - sm(ujT + 6) knu\(x)] . (27) 

If the general solution can be written as a sum over such modes is a subtle question, but here we ask mainly if there 
are any growing modes with A > 0. Furthermore, if we discretize u[x,t) on a finite grid in x (but not in r), the 
general solution of the discretized field equations is then clearly given by the sum over a finite number of modes, each 
of which is exactly exponential in r. 

In defining the variables u that we evolve numerically, we begin from the all the variables with an overbar defined 
in p], and their dot and dash derivatives where necessary. On a self-similar background, it is useful to further rescalc 
these variables by powers of e T so that the background plus perturbations is still self-similar if and only if the rescaled 
variables are independent of r. If this is done, the resulting equations do not contain explicit powers of e T . The 
spacetime perturbations grow in a physical sense towards the singularity if and only if the variables u grow with r. 
Rescaled barred variables will be denoted by a circle, their rescaled dot derivatives by a tilde, and their rescaled prime 
derivatives by a hat. 

In the following, we shall use "degree of freedom" to denote a variable that can be freely specified as a function of a 
radius at the initial moment of time. In this count, a wave equation has two degrees of freedom, for each value I and 
771. In the I > 2 perturbations there are eight physical degrees of freedom, corresponding to wave equations for the two 
polarizations of gravitational waves, the three components of the Euler equation, and the continuity equation. On a 
spherically symmetric background, these generic eight degrees of freedom split into three axial and five polar degrees 
of freedom. The number of first-order variables is larger (four and seven, respectively) because in the first-order form 
of a wave equation for a variable <p, <f> itself and <j> >x are separate variables that are linked by a (trivial) constraint 
equation. 

B. Numerical methods 

In investigating the non-spherical perturbations numerically, we have to treat each value of k and I separately. (The 
equations do not depend on 777.) In practice, we work with a finite sample of values of k in the range < K < 1, and 
with I < 5. With increasing I, numerical difficulties at the center of spherical symmetry become more pronounced, 
limiting the range of I we can investigate. Fortunately the range I < 5 is sufficient to see a trend, as we shall 
demonstrate in plots. 

We are looking for mode solutions of the form (Uw, and in particular for the dominant mode, the one with the 
lowest value of A = ReA. This objective allows a number of very different numerical approaches. Three different ones 
come to mind, and we describe them first, and then summarize our experience with the last two of them. More details 
are given in the appendixes. 

1) Making the ansatz (|2q), and imposing suitable regularity conditions, one obtains a boundary value problem for 
an eigenvector U\(x) and eigenvalue A. For the I > 2 axial perturbations, regularity conditions have to be imposed at 
the center and at the light cone. For the I — 1 polar perturbations regularity conditions have to be imposed at the 
center and the soundcone, and for the I > 2 polar perturbations regularity conditions are required at the center and 
both the light cone and the soundcone (which is in the interior of the numerical domain). Modes can then be found 
in two ways: 



a) From an initial guess for A and u\(x), one can find the correct values by shooting or relaxation. In practice, 
the initial guess has to be quite good, and finding one solution does not exclude the possibility that there is another 
solution with larger A. 

b) For given A, the boundary (and possibly midpoint) conditions can be solved in terms of n free parameters. If the 
boundary value problem is well-posed, the shooting procedure must match up n variables. The mismatch in these n 
variables is a linear function of the n free parameters, and is therefore described by an n x n matrix A that depends 
analytically on A. If A(X) has a kernel, a solution u\ can be found for this A. One therefore looks for zeros of the 
complex-analytic function det A(X). This can be done by contour integrals |12(| . 

2) One can also use the equations to evolve generic initial data u(x) in r. At late times the solution will be dominated 
by the dominant mode, and one can read off A from its time-dependence. One can also subtract the dominant modes 
one after another, in the Gram-Schmidt process , in order to find subdominant modes. This is known as the Lyapunov 
method ||. 

3) In a third approach, the evolution equations are finite differenced in x but not in r. The resulting finite 
difference-differential equations can be used in two ways: 

a) With M degrees of freedom on N grid points in x, the map T : m(xj) — ► Ui. T (xj) with i = 1...M and 
j = 1 . . . N is a square matrix of size (MN) 2 . Its eigenvalues with the largest real parts should be an approximation 
to the continuum eigenvalues A with largest real part. (The lower eigenvalues will depend on the discretization scheme 
and are not expected correspond to continuum modes.) 

b) Alternatively, we can use a standard ODE integration scheme to discretize in time. Such a numerical method 
is called "semi-discrete" because with sufficiently small step size At it is effectively discrete only in x. The map 
Ta : u(0,x) — > u(A, x) for a finite interval A is again an (MN) 2 matrix. The few eigenvalues with largest modulus 
should now be approximations to the largest of the numbers e . 

We first implemented the Lyapunov method, method 2). It is the simplest method for obtaining the dominant 
mode. This method worked well for the k = 1/3 fluid S, and also for the nonspherical perturbations of the scalar 
field critical solution p0| . However, in the polar perturbations for some values of k and I, the dominant mode is a 
numerical artifact (an instability) in all finite differencing schemes that we have tried. (The nature of the instabilities 
will be discussed below.) Then one needs to look for sub-dominant modes of the finite difference equations in order 
to find the dominant physical mode. We have found that the Lyapunov method is extremely inefficient for finding 
subdominant modes. 

We then implemented both method 3a) and 3b). In method 3b) we used first order, second order and fourth 
order Runge-Kutta integration (RK1, RK2, RK4), and an implicit second-order scheme (iterated Crank-Nicholson, 
or ICN). For a small enough time steps, the differences between these methods are negligible, and we effectively reach 
the continuum limit in r. Furthermore, in this limit, the modes and eigenvalues produced by methods 3a) and 3b) 
agree up to roundoff error. Therefore, there is no advantage in method 3b) over 3a), but a higher computational cost. 

Within method 3), many ways of finite differencing in x are possible. We have used four finite differencing schemes 
in x. The differences between them remain important at all feasible values of Ax. Two of these schemes are upwind 
schemes that explicitly use the eigenvalues of the matrix A. Both are the linearized version of Godunov schemes. 
GDI, the scheme used before in ||], is first-order accurate. GD2 is a second-order accurate version. The other two 
schemes used centered differences and are second-order accurate. CD3 uses the obvious centered differences. CD4 
uses a well-known trick to deal with terms of the form (p iX + (2Z/x)</>, which can give rise to numerical instabilities 
near the center x = 0. (The 3 and 4 are just consecutive labels.) All four schemes are defined in Appendix |j. 

Some of the numerical instabilities that we see are familiar: problems at the center, in particular for high I, and 
grid modes in centered differencing. Another kind of instability was harder to understand. The continuum equations 
admit mode solutions (pq) in which u jX is discontinuous at a characteristic that is also a line of constant x, that is, at 
the light cone and/or the sound cone of the singularity. The evolution equations are hyperbolic, and these solutions 
are called weak solutions. They are discussed in more detail in Appendix pi While they are valid as a generalized 
type of solution, these modes would not arise in a collapse situation, and so we need to exclude them. Unfortunately, 
for certain values of k and I, they dominate the top smooth physical mode. 

The weak modes were not seen in the investigation pQ] of the perturbations of the scalar field critical solution, 
because there they could be discontinuous only at the light cone of the singularity, but the numerical domain was 
truncated precisely there. The same can be done for the axial perturbations of the perfect fluid critical solution, 
because there are no axial sound waves. Similarly, the numerical domain can be truncated at the sound cone for 
the spherical and I = 1 polar perturbations because they do not comprise gravitational waves. The I > 2 polar 
perturbations, however, contain coupled sound and gravitational waves. Therefore the numerical domain cannot be 
truncated at smaller x than at the light cone, and this leaves weak modes at the sound cone. 

The finite differencing schemes we use are not designed to represent weak solutions correctly, but they do of course 
have a counterpart in the modes of the finite difference system. Sometimes the numerical counterpart resembles the 
continuum mode (in particular in GDI and GD2), but sometimes it cannot be distinguished by inspection from a 



smooth mode (in particular in CD3 and CD4). The only certain criterium is convergence. This makes it crucial that 
we have more than one finite differencing scheme, so that we can carry out independent residual evaluations, as well 
as simple convergence tests. 

We have loosely referred to the numerical artifacts as instabilities. However, the usual concept of the stability of a 
numerical method is at most exponential growth of numerical solutions. But here we are using semi-discrete methods 
(discrete only in a;) on a system of linear equations with r-independent coefficients. Therefore all solutions, both 
physical and artificial, depend exactly exponentially on r. Exponential growth, or its absence, can therefore not be 
used to distinguish between physical and unphysical solutions. An unphysical mode may grow either more or less 
rapidly than a physical mode. The only certain way of distinguishing them is by convergence. (This will also rule out 
the weak modes, because the numerical methods were not designed to handle them, and will therefore fail to converge 
on them.) 

We have carried out two kinds of convergence test. One test is to check that the discretized eigenvectors u\{xj) 
and corresponding eigenvalues A converge with increasing resolution in x at the expected rate (to first order or second 
order) for each scheme, but also that all four schemes converge to the same solution. The other test is independent 
residual evaluation. Writing the system of continuum equations formally as u.. T (x) = Lu(x), where L is a linear 
derivative operator, let L\ and L<i be two different finite difference approximations to L, and let (ui, Xi) and («2, A2) 
be modes (eigenvectors and eigenvalues) of L\ and Li. If these are approximations to a continuum mode, the norms 
|AiUi — L,2U\\ and IA2U2 — L\ui\ should converge to zero with increasing resolution. Convergence should be to second 
order in resolution if both methods are second-order accurate, and to first order in resolution if one or both methods 
are only first-order accurate. For the norm || we choose the Z 2 -norm divided by the number TV of grid points, which 
is an approximation to the L 2 norm. 

Imposing the constraints ( pa ) poses no numerical difficulty. During evolution, in methods 2) or 3b), one has the 
choice of either imposing the constraints only on the initial data, from time to time, or at each time step. We find that 
this hardly affects the results. Some care has to be taken when calculating the map T in the presence of constraints. 
This is discussed in Appendix [(J. 

C. Spherical perturbations 

The spherical perturbations have already been investigated by several authors [pH7|, an d we use them here as a 
test of our methods, and also to make sure that we are investigating the same background solution as these authors. 
Gauge-invariant perturbations do not exist for I — 0. Equations of motion for the spherical perturbations are most 
easily obtained by linearizing the field equations in spherical symmetry in the polar-radial gauge that is also used for 
the background. The perturbations 6 (hip) and SV obey evolution equations in r and x, while the perturbations 8a 
and Sg obey constraint equations in x only. We do not give the detailed equations here. 

We have used the all four codes. In GDI, GD2 and CD3 the top physical (growing) mode shows up as the top 
numerical mode. In CD4 the growing physical mode is also present but must be identified by hand because it is 
dominated by grid modes. 

Fig. g demonstrates the convergence of GDI, GD2 and CD3 with increasing resolution towards a common value 
of X(k) over the entire range of k. GDI converges approximately to first order, and GD2 and CD3 approximately 
to second order, as expected. The finite differencing error is approximately 10~ 3 in both second-order schemes at 
Ax = 1/320 for intermediate values of k. It rises towards the high and the low end of the k range. 

Figs. I and Fig. | compare our \(k) at Ax = 1/320 (GDI, GD2 and CD3) with that obtained by Maison f§. The 
results agree quite well for all k. Nevertheless, the figures show that there is a systematic difference to Maison's results 
that is generally much larger than our finite difference error. It grows in all three codes as K — ► 1. (The exception 
is in GDI at low k where the finite difference error becomes dominant.) We suspect that the systematic error is in 
our code, rather than Maison's code, and specifically in the background code: as we have discussed above, it is not 
well-behaved at the light cone as k — ► 1, and at the center as k — > 0. 
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FIG. 2. This plot demonstrates that A(k) in GDI, GD2 and CD3 converges with increasing to a common value and at the 
expected order. The resolutions are Aa; = 1/10, 1/20, ...1/320. The three plots, from top to bottom, show the error in GDI, 
GD2 and CD3. A reference value for A(fc) was obtained by Richardson interpolation on CD3 at the highest three resolutions. 
This reference value was then subtracted from the result of all three codes at all resolutions to obtain a measure of error. To 
demonstrate power-law convergence, this error was divided by a factor of 4 for each factor of 2 in Aa: (for CD3 and GD2), 
or by a factor of 2 (for GDI). The error at the highest resolution was not rescaled. Note the rise of the error at low k in all 
three methods. Although the different graphs in each plot do not lie on top of each other perfectly, they are similar, indicating 
approximate power-law scaling at the expected order. 
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FIG. 3. A(k) from Maison, and our codes GDI, GD2 and CD3 at Aa; = 1/320. (The four lines are not resolved in this plot. 



11 



0.04 



0.02 




-0.02 

0.2 0.4 0.6 0.8 1 

FIG. 4. Difference in A(fc) between Maison and our codes at Ax = 1/320: circles are GDI, squares GD2 and diamonds 

CD3. The fact that GD2 and CD3 differ from Maison by approximately the same value indicates that this difference is mainly 

systematic error, rather than finite differencing error. The fact that GD has a different deviation from Maison is explained by 

its larger finite differencing error, which increases at low k - see top plot of Fig. H. 



D. Axial I > 2 perturbations 

For I > 2 the gauge-invariant velocity perturbation [3 obeys an autonomous transport equation, just as for 1 = 1. 
For I > 2, there is also a gauge-invariant metric perturbation II, which obeys a wave equation with (3 as a source [B. 
The evolution of j3 is still autonomous. This means that the A-spectrum of the coupled II and f3 system is the sum of 
two parts: modes in which (3 vanishes, so that they are solutions of the IT equation without a source term, and modes 
that are driven by (3, so that their value of A is set by the evolution of [3. 

As shown above in section [II, the spectrum of of /3-modes can be calculated analytically for all I including I > 2. 



As a check we have implemented the (3 equation on its own, and the numerical results agree with the analytical ones, 
showing the top two of the analytically calculated modes. That is as much as one can expect, because the n-th mode 
behaves at the center as x 2n , and no finite differencing scheme is designed to represent such behavior correctly. In 
fact, unphysical modes that behave as x n with n odd also show up. They are unphysical because the corresponding 
velocity perturbation is not regular at the center, but they are valid solutions of the equation. 

The equations for II were derived in the fluid frame in || , but they look simpler in the frame of constant r observers 
(radial frame). These are just different choices of first-order variables for the same wave equation. As a test, we have 
implemented the equations in both frames, and the results converge as expected. In the fluid frame we have also 
implemented the coupling to f3. As expected this simply adds extra modes driven by (3 modes to the spectrum. As 
we have analytical results for the /3-modes (they all decay), we only need the free II modes. 

For simplicity, we give here only the source-free equations in the fluid frame. These equations are quite similar to 
the toy model wave equation of Appendix [D| The variables are 



n 



-i 



T n, 



n 



e-C'+^cT 1 !!*, 



n 



-C+^a- 1 !!, 



(28) 



II and II can be specified freely on the initial surface, while II is constrained by 

o 

II = asIT -r. 



(29) 



The II equation without source, in the radial frame, is then equivalent to 
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sg ) sg 



c 
n T = Ciir x + A x n x + — n - 

sg 

n r = -n - sxaii - ifi. 
g 



c a 



n-(i + 2)- 



i + 1 + 9± \ n . 

s 9. 



1 

s 2 ag 


'^ + (1- 

X 


a 2 -l 

X 2 


o 

n, 


(30) 






(31) 
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This is of the form u :T — Au, x + Bu, where the coefficients of the 2x2 matrix A are 

A\ = -x, Bi = — , C\ = — . 
sg sg 



The matrix A has therefore the eigenvalues 



A ld 



-x± 



*g 



(33) 



(34) 



We have used GDI, GD2 and CD3. All three codes have numerical modes that are not physical. Often it is clear 
that they are related to weak modes at the light cone, and they can be discarded by inspection. Weak modes are 
particularly easy to spot in GDI, as it is the most causal code. In the other codes, the numerical counterparts of weak 
modes can appear quite smooth. In GDI we can suppress weak modes completely by rearranging the grid spacing so 
that the light cone falls exactly on a grid point, and then truncating the grid at that point. This code will be referred 
to as GD1LC. In the other codes we need to pick out the physical modes by the criterion of convergence both with 
resolution for a single code and between different codes. Fig. || demonstrates convergence of A for the case K = 1/3 
and I = 2. The best values (either from GD2 or CD3) for the Lyapunov exponents are plotted as a function of k and 
I in Fig. [| The I — 2 leading mode is unstable for 0.58 < k < 0.87 (see also Fig. Q). 

Convergence tests fail to identify any mode as physical at the highest available resolution for k — 0.05 at I = 4 
and k — 0.1 at I = 5. However, all modes of all three codes decay for these values of k and I, so that we are fairly 
certain that there are no physical growing modes for these values of k and I. Furthermore, for all k, ReA decreases 
with increasing / (see Fig. Kj). Therefore we conclude that all Z > 3 modes are stable. 




2 3 4 

FIG. 5. Convergence of A for the top physical 1 = 2 axial (II only) mode with resolution, at n — 1/3. The upper graph is 
ReA, and the lower graph is ImA. From left to right iV = 10, 20, . . . 160. Circles are GDI, squares are GD1LC, diamonds GD2 
and triangles up CD3. 
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FIG. 6. Best value (using CD4) for the Lyapunov exponent A of II, against k and /. The upper graph shows ReA against k, 

and the lower graph ImA. Circles denote 1 = 2, squares / = 3, diamonds I — 4 and triangles 1 = 5. Points are linked by straight 

lines. The points k, — 0.05, 0.1, I — 5 and k = 0.05, 1=4 are missing because A could not be computed. The curves for ImA end 

at small k where the modes become real. 
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FIG. 7. Detail from Fig. | showing ReA for I = 2 and I 
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2 leading mode is unstable for 0.58 < k < 0.87. 



E. Polar I — 1 perturbations 



In this section, we discuss the equations of motion for the I = 1 (dipolc) polar perturbations. There are no dipole 
gravitational waves, and therefore the gravitational field has no degrees of freedom independently of the matter. 
There are three matter degrees of freedom. They are an azimuthal fluid velocity perturbation a, a radial fluid velocity 
perturbation 7, and a density perturbation u). (For the equation of state we use here, the pressure and density 
perturbations are related by Sp — n8p). The metric perturbations are obtained from the matter perturbations by 
constraints. For a general discussion we refer the reader to 0. Here we only carry out the reduction of the equations 
given there to a self-similar background solution. 

We use the matter variables 



a 



r a, 



7=— (l + «)7, uj = lu. (35) 

The matter perturbations are regular at the center if a and 7 are 0(1) and even in powers of x, and to is 0(x) and 
odd. The leading orders of a and 7 are additionally constrained as (1 + n)a + 7 = 0(x 2 ). 

The perturbation variables are not completely gauge- invariant for 1 = 1. We fix most of the gauge freedom by 
setting the metric perturbation k = 0. In this (partially fixed) gauge the remaining metric perturbations ip, x & n d r\ 
are determined by constraints. We introduce suitably rescaled metric perturbation variables 
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a 



X = e r Xi ip = e r tp, rj = e r rj. (do) 

o 

The metric perturbations are regular at the center if \, rj and ip are all 0(1) and even in x. The perturbed spacetime 
is CSS if and only if all perturbations are independent of r. 

From Eqs. (99,100) of ||, the equations of motion of the matter perturbations are 

-S-y (37) 

-S u (38) 

-r~ l S a - Ua, (39) 

where the expressions 5 T , 5 W and S a are given by Eqs. (A9-11) of M. The first two equations constitute a first-order 
form of a wave equation whose characteristics have velocity ^/k with respect to the background fluids: they describe 
sound waves. The third equation transports a along the fluid. Solving the three matter equations for the r-derivatives, 
on a continuously self-similar background spacetime, one obtains the following system of evolution equations: 

7,r = A^. x + B K & >X + 5(7), (40) 

<S,t = C*7,» + 4Ax + £(<&)> (41) 

a, T = Fa, x + S{a), (42) 



The coefficients of the 3x3 matrix A in ( 24 ) are 



(1 - k)V k(1-V 2 ) _ (1-V 2 ) V_ 

A -~- x - {1-kV 2 W Bk ~ (l-^V^sg' K ~ {l-nV 2 W F_ ^ _ ^' (43) 

Its eigenvalues, the characteristic speeds, are 

(l- K )V y/ZQ. - V 2 ) V 

\ K ±=-X- T ± , \ = -x . (44) 

(1 — nV z )sg (1 — nV z )sg sg 

The constraint equations for the metric perturbations are also given in H. The derivative operator D becomes a 
simple partial derivative in coordinates x and r, namely rD = xd/dx. We obtain a system of three ODEs in x, from 
now on referred to as the constraints: 

o 

xu :X = Mu + s, u=(x 1 r),ip). (45) 

M and s are given in Appendix pL 

We now return to the issue of the gauge freedom that is left after one has set k = 0. The change of k under an 

OOO O j—, 

arbitrary gauge transformation parameterized by the variable £ is k — > k + L\^ (see |9j for details.) L\ is a differential 
operator (roughly speaking a second x-derivative), whose kernel is parameterized by one free function of r. The 
remaining gauge freedom can therefore be fixed by imposing a gauge condition worth one function of r. There are 
two natural choices. In "77 gauge" we set rj — 0(x 2 ) at the center. In "a gauge" we set a = 0(x 2 ), and therefore 
also 7 = 0(x 2 ). Numerically it is not easy to enforce that any variable behaves as 0(x 2 ) at the center. Therefore we 
impose the two sub-gauges using variables rescaled by suitable powers of x. This is discussed Appendix pi 

Although by function counting either a gauge or 77 gauge should fix the residual gauge freedom in k = 0, a single 

o 

gauge mode does in fact survive in each gauge. The change of rj under a gauge transformation is rj — > 77 + L2<£- Any 

OO o 

residual £ must now obey both Li£ = 0, from the gauge condition k — 0, and L-£ = 0{x 2 ), from the gauge condition 
r) = 0(x 2 ), where L2 is another differential operator (roughly speaking the wave operator). A careful calculation 

o 

shows that the two joint equations still have one non-zero solution, which is £0 = e AT /(*) f° r 

(46) 

That means that in the 77 = 0(x 2 ) sub-gauge of the k = gauge we have a single complex conjugate pair of gauge 
modes left. 

o 

The situation is similar in a gauge: The change of a under a gauge transformation is a — > a + ^3^ for some 

o o 

operator L3, roughly speaking a first r-derivative. The joint kernel of L\^ = and L^^ = 0(x 2 ) is not empty either 
but contains a single real gauge mode with 
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A = 



l + 3fc 
3(1 + A)' 



(47) 



Bothgauge modes are in fact found numerically, and give us a check on the numerical precision for each value of k. 
Fig. |8] demonstrates this for both gauges. 

We have used GDI, GD2 and CD4. For k > 0.15, the numerical results are straightforward. Both in a gauge and in 
r\ gauge, we find the expected gauge mode. We find that a gauge, in the specialized variables described in Appendix 
pi, has the best numerical behavior at the center. As an estimate of numerical error in the code, Fig. shows the 
exact and numerical value of A for the gauge mode in a gauge. Fig. Uu shows the best value of A in this system as a 
function of n. 

For k < 0.15, our code does not work well. Fig. g shows that with decreasing k, the error in locating the gauge 
mode increases rapidly. For k = 0.1 and k = 0.05, all three codes fail to find the real gauge mode. Nevertheless, 
GDI and GD2 have precisely one growing mode, and CD4 has precisely one growing mode that is not clearly a grid 
mode. All these modes are complex, but they obey |Reu(a:)| 3> |Imu(x)| for all x, by a factor of ~ 100: the mode is 
almost real, modulated by a coswr factor. It is also clear that there is a finite differencing problem at the center that 
affects all three codes: the functions u(x) are not well-behaved even or odd functions. The origin of this instability 
is uncertain. A plausible explanation is that the numerical background solution itself is not sufficiently well-behaved 
at the center for small k. Some of the coefficients in the I = 1 (and also I > 2) polar perturbation equations are 
very large and sharply peaked at the center. Furthermore, some of these coefficients must be obtained as numerical 
derivatives of the background solution. These numerical derivatives are not smooth at the center, and we had to 
artificially extrapolate them to the center in order to smoothe them. 

Because all three codes agree on a single growing mode that is almost real in the sense just defined, we identify 
this mode with the gauge mode, even though the agreement in A is poor. This leaves no other growing modes in any 
of the codes for any k (except for obvious grid modes in CD4 that can be ruled out by inspection.) We therefore 
conjecture that there are no growing physical modes, even though we cannot identify the top physical mode. 
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FIG. 8. Convergence of A for the single gauge mode in the I = 1 polar perturbations, for k = 1/3. From left to right 
N = 10, 20, . . . 160. Circles are GDI, squares are GD2, and diamonds CD4. The top two graphs were obtained in r/ gauge, 
using the variables of Appendix A, and show ReA and ImA for the gauge mode. The value of A for this mode computed from 
the background solution and Eq. (|4£]) is A ~ 1/2 + 2.384i for k — 1/3. The bottom graph was obtained in a gauge, using the 
variables described in appendix B. The exact value of A for this gauge mode is 1/2 for k — 1/3. 
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FIG. 9. A for the single gauge mode in the I — 1 polar perturbations in a gauge as a function of k. Circles are measured 

values, connected by straight lines. The smooth curve is the exact value. The points at k = 0.1 and k — 0.05 represent the real 

part of the only growing mode present, which is complex. 
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FIG. 10. Best value (obtained with GD2) for the physical Lyapunov exponent A of the I = 1 polar perturbations. The 

upper graph shows ReA against k, and the lower graph ImA. Points are connected by straight lines. The measured values end 

at k — 0.15 because for k = 0.1 and k — 0.05 the top mode could not be identified. The graph of ImA ends because A is real 

for k < 0.5. Clearly a real mode and a complex mode pair have crossed between k — 0.5 and k = 0.6. 



F. Polar / > 2 perturbations 

As discussed in detail in ||, the / > 2 polar perturbations of a spherical perfect fluid admit a free evolution 
formulation. The dynamical variables that can be specified freely are the metric perturbations x, k and ijj, and 
the time derivatives x an d k. \ an d k obey wave equations, and ip a transport equation, all of which are coupled. 
There are three gauge-invariant matter perturbations, a is an azimuthal velocity perturbation, 7 is a radial velocity 
perturbation, and w is a pressure and density perturbation. These matter perturbation variables are completely 
determined as algebraic expressions in the seven first-order metric perturbation variables, as discussed in detail in [g| . 
For the purpose of finding the late-time behavior of generic perturbations, we work with the metric perturbations 
alone. 

Solving for the r-derivative of the seven variables, we have find that the matrix A in (p3) takes the following form: 
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A = 



( A x B x D \ 

1 d A x vD l 

^k B K E 

C K A K vE 

V F J 



X 

k 
k 



(48) 



The remaining two variables \ an d k are evolved using only source terms. The coefficients of the matrix A have 
already been given in Eqs. (|3|) and (E3h, except for 



D 



9 

s z gx 



E = - 



2kxU(1 - V 2 ) 



(49) 



The five eigenvalues of A are Ai±, X K ± and Ao, which have already been given. (D and E do not influence the 
eigenvalues.) Note that s has been chosen so that A K + > for < x < 1. The point x = 1 where it changes sign 
is the sound cone. Similarly, the point x = x\ c > 1 where Ai+ changes sign is the light cone. Clearly X\ c is 1 for 
n = 1 and diverges as k — » 0. Note that the (x, \) sub-matrix of A is the same as the entire matrix A for the axial II 
perturbations, as both pairs of variables obey a wave equation at the speed of light. Similarly, the (fc, k) sub-matrix 
is equal to the entire matrix for the I = 1 polar perturbations, which obey a wave equation at the speed of sound. 
Therefore A K is equal to Ax for k = 1, when the speed of sound is the speed of light. 

We have used the three codes GDI, GD2 and CD3. There is good agreement and convergence between the three 
codes for 0.2 < k < 0.8. Fig. |ll] demonstrates convergence of A for k = 1/3 and / = 2. Fig. [lj shows the best values 
of A for I — 2 ... 5. We find that all modes for all I decay for n < 0.49. For k > 0.49, there is an unstable I = 2 mode. 
At still higher k, there is probably more than one unstable mode, but it is difficult to identify subdominant modes 
with certainty. 




1 2 

FIG. 11. Convergence of A for the top I = 2 polar mode with resolution at k = 1/3. From left to right N — 20, 40, . 
Circles are GDI, squares are GD2, and diamonds CD3. The upper graph is ReA, and the lower graph is ImA. 
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0.2 0.4 0.6 0.8 1 

FIG. 12. A for the top I — 2. ..5 polar modes as a function of k. The upper graph is ReA, and the lower graph is ImA. Circles 

are I = 2, squares I = 3, diamonds I — 4 and triangles I = 5. Note that ReA decreases with I, while ImA increases. The curves 

for / = 3, 4, 5 do not extend to the lowest and highest k because of large numerical error. 



0.5 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 

FIG. 13. Detail of Fig. [l2l ReA is positive for I = 2, k > 0.49. ReA is negative for all k for I > 2. 



For k > 0.8, and k < 0.2, the numerical precision decreases quickly. In this region the numerical accuracy is poor 
for I = 2, and for I = 3, 4, 5 even the top physical mode cannot be identified clearly. The origin of these problems is 
different for high and low k. At high k, all modes in all three codes show discontinuities at the sound cone. These 
are not due to the weak modes discussed above, but are due to unsmooth behavior in the numerical background 
solution, through certain coefficients of the perturbation equations that must be obtained as numerical derivatives of 
the background solution. At low k, the origin of the low precision is not so clear. It may be that the background 
coefficients become increasingly large and sharply peaked in x near the center as k ~> 0. With the sound cone at 
x = 1, the value of x at the light cone also diverges as 1/vk as the sound speed goes to zero. (Recall that x is defined 
so that the sound cone is at x = 1.) In a second source of error, some of the perturbation coefficients, as determined 
from the numerical background solution, must be smoothed at the center for small k, and this introduces additional 
error. (This is the same problem that affects the / = 1 polar modes.) 

In contrast to the inadequacy of the numerics at high and low n, the numerical precision is good in a large 
neighborhood of the value k ~ 0.49. It is therefore certain that the I = 2 polar perturbations are stable for k < 0.49 
and unstable for n > 0.49. As evidence that the instability is physical, Fig. [14] demonstrates the convergence of 
independent residual evaluations between all three codes for I = 2 and n = 0.6. From convergence, we estimate A (for 
this k and I) as A = (0.112 ± 0.003) + (1.968 ± 0.005)«. The real part of A is therefore much larger than the finite 
differencing error. 
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12 3 4 

FIG. 14. Log-log plot of independent residual evaluation between three different codes for the dominant polar I — 2 

mode for n = 0.6. Evaluating the mode produced by GD2 with CD3 shows exact second-order convergence (triangles up, 

bottom curve). The opposite operation, checking on CD3 with GD2, shows second-order convergence that breaks down at high 

resolution (triangles down). Checking on CD3 and GD2 with GDI shows the expected first-order convergence (triangles left 

and diamonds). The opposite check shows first-order convergence that breaks down at low resolution (squares and circles, top 

curves). The resolutions shown are Ax — 1/10 . . . 1/160. 

V. CONCLUSIONS 

The threshold of gravitational collapse of a perfect fluid with p = kp in spherical symmetry shows "type II" critical 
phenomena (including mass scaling with a critical exponent) for the entire range < k < 1 Mn- Is this true also when 
one relaxes the assumption of spherical symmetry? We have addressed this question by studying the nonspherical 
perturbations of the critical solutions in spherical symmetry. We have extended the results of H from k = 1/3 in 
the equation of state p — kp to the entire physical range < k < 1. Because of the numerical difficulties, we have 
complemented the single first-order scheme used in M with three different second-order finite differencing schemes. We 
have complemented convergence tests by independent residual evaluation. This allows us to identify "weak modes" , 
numerical artifacts related to weak solutions of the perturbation equations. Finally, we have corrected errors in jq| 
concerning the equations of motion for the I = 1 axial and polar perturbations. 

After verifying that the CSS solution has exactly one growing spherical perturbation for all k, we have found that 
all nonspherical perturbations decay for all k with the following exceptions: 

1. For k, < 1/9, there is precisely one growing I = 1 axial mode. This result was obtained analytically, even though 
the background solution is known only numerically. Because it is I = 1, this mode is three-fold degenerate. 

2. For k > 0.49, there is a growing I = 2 polar mode. We cannot rule out that at larger k, there are several growing 
modes. Because it is / = 2, this mode is five-fold degenerate. 

3. For 0.58 < k < 0.87, there is also a growing I = 2 axial mode. 

The numerical evidence for this is not as good as one would like. Because of numerical error, we cannot measure 
A of the top physical mode of the I = 1 polar perturbations for k < 0.15, nor of the / > 3 polar perturbations for 
k < 0.2 and k > 0.8. However, in all these cases where physical modes and numerical modes cannot be distinguished 
clearly, all modes, including the physical modes, do in fact decay. Therefore we still argue that all physical modes 
decay. One cause of the numerical difficulty is the separate existence of light cones and sound cones, which gives rise 
to numerical artifacts related to weak solutions of the continuum equations. A second cause is that instead of a single 
smooth background solution we are dealing with a one-parameter family of such solutions, which is ill-behaved at 
both ends k = and k = 1. A third cause is that several coefficients required in the perturbation equations are first 
and second derivatives of the background fields, which are not perfectly smooth at the center and the light cone. 

What is the significance of our results? The I = 1 axial perturbations are naturally associated with infinitesimal 
(differential) rotation. The presence of a growing rotation mode at low k (< 1/9) is not surprising, as one would expect 
a rotating dust configuration to be torn apart by centrifugal forces. For a sufficiently stiff fluid this intuition obviously 
fails. The significance of the I — 1 axial perturbations is that they can survive into the final black hole formed in 
collapse (and turn a Schwarzschild black hole into a Kerr black hole) , while all other non-spherical perturbations must 
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be radiated during collapse. Their instability is expected to give rise to interesting new phenomena in critical collapse 
that will be explored elsewhere |L3] . 

The existence of a small number of unstable polar modes at high k (> 0.49) is more puzzling. At face value it 
suggests that the spherically symmetric CSS solution is a critical solution only when restricted to exact spherical 
symmetry. Its place may be taken by another critical solution with less symmetry or there may not be a critical 
solution, and hence no universality, at the black hole threshold. 
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APPENDIX A: BACKGROUND EQUATIONS 

The relation between the fluid frame derivatives and the partial derivatives in CSS coordinates is 



/ = 



/' = 



eg 



ay/1 - V 2 
e T g 



■I' + [ X+ Jj) ' 



ay/1 ~ V 2 
rDf = xf yX . 



Vf, T 4 ( Vx + — ) f,a 



(Al) 

(A2) 
(A3) 



The perturbation equations of |9| allow for a 2-parameter equation of state p = p(p,s), where s is the entropy per 
particle. When restricting to our simple barotropic equation of state, we set both the fluid entropy s and its gauge- 
invariant perturbation a to zero. We also set C = dp/ ds to zero, and we set the sound speed squared c 2 = dp/dp 
equal to the constant n. 

The background equations that result from the CSS ansatz are 



9 , 2s X p 
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Here (A5) and (A6) are the fluid equations of motion, and (A7) and (A4) are two of the Einstein equations. The 
terms C a and C a arise when metric derivatives in the matter equations are eliminated using the Einstein equations. 
The expression for a incorporates the regularity condition (absence of a conical singularity) a = 1 at r = 0. 
The following background quantities are required as coefficients in the perturbation equations. 



W = e~ T W = e- T n A r A = (asx)' 1 ^ - V 2 y 1/2 = 0(x~ l ), U = e- T U = e~ T - A 



u~r, A = VW = 0(1), (A13) 



p, = e T p = e T u A \ A = 0(l), D = e T v = e T n A i A = O(x), 



p=e- 2T ATTp^O{l), 



2m 



1 - r iA r A = 1 - r 2 (W 2 - U 2 ) = 1 - a~ 2 = 0(x 2 ). 



(A14) 
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In a CSS background all these quantities depend only on x, and we have indicated their behavior at the origin. 



APPENDIX B: L = 1 POLAR PERTURBATION EQUATIONS 

The source terms in Eq. ( |4C| ) are obtained from the equations of [^| in the form 
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where the terms on the right-hand side are defined in Ref . B] . Fully expanded in the variables adapted to self-similarity, 
they are 
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The coefficients in the constraint equations Eq. (|45|) are 
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(B7) 

s = a 2 ( -2p s f + 4p{l + K)W(sx)a 1 (1 - n)p(sx)u) - 4(1 + K)pU(sx) 2 a, 2p(sx)- 1 & - 8(1 + n)pUa\ . (B8) 

In u r) gauge" we set i] = 0(x 2 ) at the center. The constraint equations are then solved with the boundary conditions 

o 9 3? 

1> = 2(1 + K)p&, ?i = 0, X =-p(sxy 1 £ J + 6 --p a (B9) 
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at the center x = 0. The condition rj — is our gauge condition, and the other two conditions follow from Mu + $ = 
(M has rank 2 at the center). An alternative way of imposing r\ gauge is to introduce the new variable 

fj = (sx)- 2 rj, (BIO) 

which i s 0(1) at the center in this gauge. The evolution equations are as before, with only r\ replaced by (sx) 2 fj in 
(B4-Bq). The constraint equations become 



xu, x = Mu + s, u = (x,ry,^), 



with 



-1 + a 2 [-2 + (1 - k)(sx) 2 p] a 2 (4UW +_2pW - 2vU)(sxf 2_a 2 (vU_ - pW)(sxf 
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s = a 2 (-2p7 + 4p(l + n)Wsxa, (1 - k)p(sx)- x 5j - 4(1 + n)pUa, 2p(sxy 1 u - 8(1 + n)pUa). 
The boundary conditions at x = are 

ij) = 2p(l + n)a, fj = --p(sx)- 1 uj + — pa, x = -p(sxy 1 u;+ — pa, 

from Mu + s = (M has rank 3 at the center). 

In "a gauge" we set a = 0(x 2 ) (and therefore also 7 = 0(x 2 )). The constraints are then solved with the boundary 
conditions 
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y 9(i + K ) 3 , 

Here the boundary condition on r\ comes from consistency with the evolution equation a iT = 0(x 2 ). An alternative 
way of imposing a gauge is to introduce new matter variables 

lu = (sx)~ 1 £}, 7 = k~ 1 (sx)~ 1 j, a= (sx)~ l a = sxa. (B16) 

The evolution equations in these variables are 

w )T = A K u!, x +B K j tX + S(u>), (B17) 

7,r = C K tD, x + A K % + Stf), (B18) 

a, T = Fa^ x + S(a). (B19) 

Note that the coefficients of the x-derivatives are unchanged, but that uj has taken the place of 7 and 7 has taken the 
place of ui. Accordingly, 7 is now odd and 0(x) while uj is even and 0(\). The source terms are 
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The variable a, which is odd and O(x), is much better behaved than the more obvious definition a = (sx)~ 2 a. (The 
combination fj + k/(1 + k)lu is 0(x 2 ) at the center. Numerically, it is easier to divide this by sx in S(a) to obtain 
an 0(x) term than to divided by (sx) 2 .) The constraints are solved with the matrix M given above and the source 
terms 

s = a 2 {-2pK(sx)j + Ap(l + K )W(sx) 2 a, (1 - k)p(sx) 2 uj - 4(1 + K)pU(sx)' i a, 2pu> - 8(1 + n)pU(sx)a). (B23) 



The boundary conditions are (B15) with (sx) d replaced by u> 



APPENDIX C: L > 2 POLAR PERTURBATION EQUATIONS 
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In order to work with variables that are regular and O(l) at the origin for any I, we redefine them as \ 
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ij) and k = r ■ k, as in ||. In order to obtain a first-order formulation, we introduce the frame derivatives of 



the barred quantities with respect to the fluid frame. From Eqs. (87-89) of |9J, the equations then take the form 
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(k'y- (%)' = v%- pk' = e {l+2 ^ T S A , (C4) 

W - -r-P+^Sv, - (I + W ee e« + V T S 5 , (C5) 

where the source terms S x , S k and S^ are given in Eqs. (A1-A3) of B|. The second and fourth equations are just 
the commutation relations between the dot and prime derivatives. Here they serve as auxiliary evolution equations 
for x' and k! . To these five equations we add the trivial evolution equations (%)' = x an d (k)' = k. 

For our particular application to a self-similar background we further rescale these first-order variables to obtain 
the final dynamical variables 
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The variables with a circle or a tilde are even and O(l) at the center, while the variables with a hat are odd and 
0(x). The perturbed spacetime remains CSS if and only if all seven perturbation variables are independent of r. 
These seven variables obey evolution equations without any c onstraints , except for the trivial ones that arise when 
one writes a wave equation in first-order form, given below in (C15 C16). Applying the same rescaling to the source 

o 

terms S\ to S$, we obtain, in the notation Bu = S(u), the source terms in the evolution equations for x to ip. To 

o 

these we add two evolution equations y and k, which have only source terms, but no ^-derivatives. We obtain them 
by solving (Al A2) for / T in terms of / and /'. Putting all seven source terms together, we have 
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Finally, solving ( |Al| , |A2| ) for f_ x in terms of / and /', we obtain two constraint equations, that is, equations that do 
not contain r-derivatives. They are 



X,x = 

O 

fc.X = 



as 



ix-vx), 

(k-Vk). 



Vi -v 2 

The intermediate source terms Si are 
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We have used the background equations of motion, and have introduced the new coefficient e~ 2r (p' — u), which is 
again independent of r on a CSS background. Note that all coefficients in the evolution equations are explicitly 
regular when one takes into account that 77 — fx = 0(x 2 ), p — (m/r)/(sx) 2 = 0(x 2 ) and e~ 2T (p' — v) = 0(x). The 
two terms 2(1 + l)Wx and 2k(Z + l)Wk are the equivalent of the term 2(1 + l)</>'/r of the toy model wave equation 
that we discussed above. They can be regularised in the same way. 

APPENDIX D: A TOY MODEL FOR THE EVOLUTION EQUATIONS 

We have tested various numerical methods on a toy model, the scalar wave equation in flat spacetime. This model 
is also useful as an illustration of the types of variable and the methods we use. 

Let <I> obey the free wave equation on the flat spacetime ds 2 — —dt 2 + dr 2 + r 2 dfl 2 . We make the ansatz 



$ =^2<t>im(r,t)Y lm ($,cp). 

l.rn 

We now consider a particular value of Z and m, and drop these suffixes. Then 4n m obeys 
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As in Pi we introduce the rescaled and first-order variables 
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$(£, r, 9, if) = $(i, x, y, z) is analytic in Cartesian coordinates at the origin if and only if </> is analytic in r with only 

even powers of r. In this sense <f> and <fi are even functions of r, and are generically finite and nonzero at r = 0, and 
(j)' is odd and generically 0(r) at r = 0. 
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We now go over to self-similarity coordinates x and r denned in ([I]) (with s = 1). In order to mimic the fluid 
perturbation equations, we rescale the first-order variables once again as 

l = e- nT 4>, j>=e- {n+l]T 4>, 4>=e- {n+l)T 4>'. (D4) 

The most natural choice of n for the toy model is n — 0, but for the fluid perturbations, it will be fixed by the 
requirement that the perturbed spacetime remains self-similar if the perturbation variables do not grow or decay with 
r. We therefore leave it free. We finally obtain the first-order system 

~4> T = -x4> x + x + 2(1 + 1)- - (n + 1)4>, (D5) 

x 

4>,t = 4>,x - %4>,x - (n + l)<j>, (D6) 

o o 

<p T = 4> — X(f) — n<f). (D7) 

There is also the trivial constraint 

°<P, X = $ (D8) 

that follows from the definition of the first-order variable <f>. 

o 

One significant property that the toy model has in common with the full perturbation equations is that <f> and <p 
are even and generically nonzero at x = while <f> is odd and generically 0{x). All terms except (fr/x are explicitly 
regular, and do not require special treatment at the origin. In particular, there is no term of the type 4>/x 2 left. 
Experience shows that terms of the form <p/x, where </> is an odd function of x and 0(x), can be regularised without 
giving rise to instabilities, while division by x 2 is much more troublesome. 

In other, unimportant, aspects the toy model is simpler than the full perturbation equations. As there is no fluid 
frame in this toy model, we choose the dot and prime to be frame derivatives with respect to the constant r frame. 

o 

As the spacetime is flat, the dot and prime derivatives commute. Finally, in the toy model (f> does not couple back to 
<f> and (j>, and therefore plays only a passive role. 

APPENDIX E: FINITE DIFFERENCING 

For the evolution equations (E4) at hand, the eigenvalues and eigenvectors of the matrix A can be calculated 
in closed form. The eigenvalues of A are dx/dr on characteristics of the equations: fluid world lines, radial light 
rays, and radial matter (sound wave) characteristics. The characteristics are symmetric around the line x — 0, but 
with increasing x they tip over until at sufficiently large x all eigenvalues are negative. This means that at large x 
information travels only from smaller to larger x. The reason for this is of course that while x = and lines of small 
constant x are timelike, lines of large constant x are spacelike. The "outer boundary" x = x max of our numerical 
domain < x < x max is therefore a future spacelike boundary, and so no boundary condition is required there. 

In order to use it to obtain a free boundary condition, we make the numerical method reflect the propagation of 
information, so that at x — x max all ^-derivatives are calculated using one-sided finite differences. Following jlj], we 
split A into a left and a right-moving part. Let V be the matrix of (column) eigenvectors of A. Let A be the diagonal 
matrix composed of the corresponding eigenvalues. Then A — VAV~ l . Let A + be A with zeros in the place of the 
negative eigenvalues, and let A_ be A with zeros in the place of the positive eigenvalues. Then define A± = VA±V~ 1 . 
It is clear that A = A + + A_. We now use A + with left derivatives and A_ with right derivatives. 



dr 



{A+D + u] + A-D-u] + Bu]), (EI) 



with the one-sided derivatives 



u'* . i — u" u- — u' i 

D + u] = i +1 \ D_u] = > '- 1 . (E2) 

J Ax J Ax 

Here the coefficient matrices A±, B and C are evaluated at Xj. After each time step, the constrained variables w are 
obtained at r n+1 from the u n+1 . This scheme is first-order accurate in x. (Note that we have not discretized in t 
yet.) Left differences at x = are evaluated using ghost points based on the fact that all grid functions u are either 

2G 



even or odd in x. This scheme can be thought of as the Godunov method applied to a linear equation, and we shall 
refer to it as the first-order Godunov scheme, or GDI. To obtain a second independent finite differencing scheme, we 
have also used the second-order one-sided derivatives 



2Aa; ' * ~ 2Ax 



D + u 3 ^Z > D - u 3 KTZ ■ ( t3 ) 



We shall refer to this scheme as GD2. At the outer boundary x — x max , no right derivatives are required, as there all 
information travels from the left to the right. (A + vanishes.) All variables u are either even or odd functions of x. At 
the inner boundary x = 0, left derivatives are calculated using fictitious grid points at negative x that are obtained 
as u(—x) = ±u(x). By construction, the Godunov method has the advantage that it does not require a special outer 
boundary conditions. Unexpectedly, it also has the advantage that it handles terms of the form cj>/x term at the 
center without any special treatment. It is also completely free from high-frequency grid- modes. This last property 
is less surprising when one thinks of GDI as centered differencing plus a dissipative term |14fl. 

We have also implemented a more standard finite differencing scheme based on centered differences in x, with a 
(first or second order) one-sided derivative at the outer boundary x = a; max - The center is again handled by ghost 
points, using u(—x) = ±u(x). Using naive centered differences, the code is unstable at the center because of the 
source term 2(1 + l)(f>/x in the toy model, and similar terms in the perturbation equations. A well-known remedy is 
to include this source term into the transport terms in the way suggested by the identity 



2(Z + 1) ^ _ /07 , ,d(x 21 ^ 2 
d(x 



* = (« + 3)^df (E4) 



In the toy model wave equation this procedure slows down the central instability enough so that it can be suppressed 
by numerical viscosity. We have added a centered difference expression for u, r = cu_ xx + . . ., with c of the order of 
10 -3 . However, numerical viscosity falsifies the results too much both in the toy model and in the actual problem, 



and has therefore not been used in any of the results of this paper. We shall refer to centered differencing with (E4) 
as CD3, and without as CD4. 

APPENDIX F: WEAK SOLUTIONS 

We now discuss a problem that affects the fluid and metric perturbations that we want to investigate, but that is 
already present and more easily understood in the toy model. 

The spacetime point r = t = is singled out in the self-similarity coordinates x and t, but it is not preferred on 
the flat background. A solution arising from generic C 2 initial data at t — to < has finite t and r derivatives at 
t = r = 0. The coordinates x and r, however, have been designed to "zoom in" on this spacetime point: we are 
looking at a smooth solution on ever smaller scales. Therefore, as t — > oo, a generic solution of the toy model wave 
equation should behave at large t as 

0->0(O,O)e-" T , 0->£ t (O,O)e-< n+1 > T , 0->i0, rr (O,O) e -("+ 2 ) r . (Fl) 

In deriving this fall-off we have assumed that <f> is at least twice differentiable everywhere. But as x = 1 is a 
characteristic of the wave equation, </> and its derivatives are allowed to be discontinuous there. In particular, data 
on x > 1 do not influence the solution oni< 1. Therefore, solutions exist that vanish on x < 1 for all r but not for 
x > 1. In particular, making a power-series ansatz for the region x > 1, 

4>{x,T) = e Xr [<Ma;-l) + 2 (a;-l) 2 + ...] , $(x,t) = e Xr [^(i - 1) + fa{x - l) 2 + . . .1 , (F2) 

we find a formal solution with A = I — 1 — n, which is slower than the expected falloff by I. We believe that a 
corresponding weak solution of the wave equation exists that is analytic for x > 1 (and vanishes for x < 1), but 
have not proved it. A finite difference counterpart certainly does exist, and for GDI and GD2 it falls off (or grows, 
depending on n and I) with the calculated value of A. It therefore dominates over the generic smooth solution at 
large r. We want to exclude it in studies of critical collapse, as we are interested only in the time evolution of smooth 
perturbation initial data. Numerical schemes, however, do not distinguish between smooth and unsmooth data, and 
this solution turns up in some of them, hiding the everywhere C 2 solutions we are interested in. In the future we 
shall refer to these as "weak solutions at the light cone" or simply "weak modes". We should stress that, depending 
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on the finite differencing scheme, weak solutions of the finite difference equations may or may not converge to weak 
solutions of the continuum equations, but they are always there. 

We find that both Godunov schemes when applied to the toy model wave equations develop weak solutions that are 
continuous but not differentiable at x = 1. For GDI, the finite difference mode corresponding to the weak solution 
can be analyzed easily. u tT at the first grid point with x > 1 depends only on u at that point, and on the next point 
to the left - but there all fields vanish in a weak solution. A calculation shows that u at the first grid point with x > 1 
depends exponentially on r. The exponent is either A = I — 1 — n + O(Ax), which is the continuum weak solution, or 
A = — 2/Ax + O(l), which is a finite differencing artifact that decays quickly. 

In GDI, the dominant weak is clearly not differentiable at x = 1, and this also shows up in independent residual 
evaluation as a sharp peak. If the last grid point is exactly x = 1, the numerical domain can consistently be truncated 
there, and no weak modes can arise. Because of its wider stencil GD2 is not completely causal, and so part of the 
weak mode propagates to x < 1. The numerical equivalent of the weak mode therefore appears differentiable at x = 1, 
but there is still a peak there in independent residual evaluation, although less sharp. If we attempt to truncate GD2 
at X = 1, we need to introduce an explicit boundary condition there. Doing this by using a first-order right derivative, 
or a centered derivative at x = 1 — Ax introduces spurious reflections at the boundary. Without truncation, we have 
tried updating the grid point x = 1 + Ax by interpolation from the neighboring grid points. This does not suppress 
the weak mode. Updating the three points x = 1, x = 1 + Ax and x = 1 + 2 Ax by interpolation results in a different 
spurious mode. This holds both for GDI and GD2. 

APPENDIX G: IMPOSING THE CONSTRAINTS 

Some perturbation equations contain only x-derivatives, and so are constraints rather than evolution equati ons . 



Even the free wave equation has such a constraint when it is written in first-order form. Constraints of the form ( DS ) 
are solved by integration from the center out, using the trapezoid rule, 

m = 4 + ^(> +1 + J ), (Gl) 

o 

where the starting point <fi a t the center is determined by the evolution equation. This method is second-order 

o 

accurate. We also use the exact of these finite difference equations to obtain cj) from (f> by differentiation. 

Time evolution commutes with the constraints in the continuum limit, but in the discretized equations this is only 
approximately true. In order to find the modes, we need a matrix T that acts only on a set of functions u that can 
be freely specified in the initial data and is therefore of full rank. We start with a larger matrix T' and have to 
reduce it using a numerical solution of the constraints. This reduction, however, is not unique, and although different 
reductions are equivalent in the continuum equations, they are not in the finite difference equations. 

o 

We discuss these issues in the toy model. The matrix T' acts on </>, cj) and 0, so that T" is a (3iV) 2 matrix. There are 
two natural choices for the free variables on which T acts: either <p and 4>, or the functions <j) and <p plus the number 

o ^ 

</>(0). Note that (f>(0) = by definition, so th at in either case T is a (27V) 2 matrix. 



We denote the constraint solution scheme (Gl) by / (for integration) and its numerical inverse by D (for differen- 



tiation) . In loose matrix notation we can then write the first possibility as 
and 




T* = n , n T ' 1 . (G3) 



We find that Ti works much better than T\, in having grid modes. This is not surprising giving that integration has 
a smoothing property. An alternative to T\ is to re-impose the constraints by integration after evolution, 




r 3 = . " " r p d . (G4) 
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But the matrix T3 is similar to the matrix T2, 

*-(; £)*(;?). <«) 

and so T2 and T3 have the same eigenvalues, although of course not the same eigenvectors. This is borne out 
numerically up to small rounding errors. 

In the spherical and the I = 1 polar perturbations, nontrivial constraint equations arise which are of the form 

xu x = M(x)u + s, (G6) 

where u stands for a subset of variables, and the source terms s are linear in the other variables, already known at this 
time level. The coefficient matrix M(x), an even function of x, is typically nonzero at the origin. The ODE is therefore 
singular at X = 0. We look for solutions u(x) that are even regular functions of x. These obey M(0)u(0) = s(0). 
From this boundary condition, the ODEs are solved by the implicit, second-order accurate scheme 

j/i+i = (1 - eMi+i)~ [(I + eMi) yi + e (si + s i+ i)] , e = (x i+ i - Xi)/(x i+ i +x t ). (G7) 
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